%--------------------------------------------------------------------------
% draws the eigenfunctions at time t = 0 (ASSUMES THAT h = 1)!
%--------------------------------------------------------------------------

function draw_eigenfunctions(k, p)
close all

Nx = 30;

if nargin < 2
    p = params;
end

[ev, k, M] = comp_eigs(k, p, 1);
[c, ev] = eigs(M, 1, 'SM');

fprintf('eigenvalue: %.4e\n',ev);

if (k < 1e-2)
    w = @(z) p.Ma * z.^2 .* (k^2 / 4 * (1 - z) + p.delta / p.Ma / 2 * (3 - z));
else
    A = w_coeffs(c(p.N), 1, k, p);
    w = @(z) (A(1) * z + A(2)) .* cosh(k*z) + (A(3) * z + A(4)) .* sinh(k*z);
end

x = linspace(0, 4 * pi / k, Nx);
z = linspace(0, 1, 10)';

ee = exp(1i * k * x);
for i = 1:Nx

    C(:,i) = c * ee(i);
    W(:,i) = w(z) * ee(i);
    U(:,i) = 1i / k * imag(w(z + 1i * 1e-4) / 1e-4) * ee(i);
    
end

[X, Z] = meshgrid(x, z);

plot(linspace(0,1,p.N), c);
figure;

contourf(x, linspace(0,1,p.N), real(C), 100, 'linestyle','none');
hold on;
quiver(X, Z, real(U), real(W), 'k*');
colorbar

% surf(x, linspace(0,1,p.N), real(C)); shading interp